Introduction à la modélisation statistique bayésienne

Ladislas Nalborczyk & Thierry Phénix

Planning

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 & 4: Modèle de régression linéaire

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Modèles multi-niveaux

Cours 9: Introduction à brms

Cours 10: Data Hackaton

Null Hypothesis Significance Testing (NHST)

On s'intéresse aux différences de taille entre hommes et femmes. On va mesurer 100 femmes et 100 hommes.

men <- rnorm(100, 175, 10) # 100 tailles d'hommes
women <- rnorm(100, 170, 10) # 100 tailles de femmes
t.test(men, women)

    Welch Two Sample t-test

data:  men and women
t = 2.1475, df = 197.71, p-value = 0.03297
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.2655634 6.2364837
sample estimates:
mean of x mean of y 
 173.0157  169.7646 

Null Hypothesis Significance Testing (NHST)

On va simuler des t-valeurs issues de données générées selon l'hypothèse de non-différence entre hommes et femmes.

nSims <- 1e4 # number of simulations
t <- rep(NA, nSims) # initialising an empty vector

for (i in 1:nSims) {

    men2 <- rnorm(100, 170, 10)
    women2 <- rnorm(100, 170, 10)
    t[i] <- t.test(men2, women2)$statistic

}
t <- replicate(nSims, t.test(rnorm(100, 170, 10), rnorm(100, 170, 10) )$statistic)

Null Hypothesis Significance Testing (NHST)

data.frame(t = t) %>%
    ggplot(aes(x = t) ) +
    geom_histogram() +
    theme_bw(base_size = 20)

plot of chunk unnamed-chunk-5

Null Hypothesis Significance Testing (NHST)

data.frame(x = c(-5, 5) ) %>%
    ggplot(aes(x = x) ) +
    stat_function(fun = dt, args = list(df = t.test(men, women)$parameter), size = 1.5) +
    theme_bw(base_size = 20)

plot of chunk unnamed-chunk-6

abs(qt(0.05 / 2, df = t.test(men, women)$parameter) ) # two-sided critical t-value
[1] 1.972035

Null Hypothesis Significance Testing (NHST)

alpha <- .05
abs(qt(alpha / 2, df = t.test(men, women)$parameter) ) # two-sided critical t-value
[1] 1.972035

plot of chunk unnamed-chunk-8

Null Hypothesis Significance Testing (NHST)

tobs <- t.test(men, women)$statistic # observed t-value
tobs %>% as.numeric
[1] 2.147452

plot of chunk unnamed-chunk-10

P-values

A p-value is simply a tail area (an integral), under the distribution of test statistics under the null hypothesis. It gives the probability of observing the data we observed (or more extreme data), given that the null hypothesis is true.

\[ p[\mathbf{t}(\mathbf{x}^{\text{rep}}|H_{0}) \geq t(x)] \]

t.test(men, women)$p.value
[1] 0.03297456
tvalue <- abs(t.test(men, women)$statistic)
df <- t.test(men, women)$parameter

2 * integrate(dt, tvalue, Inf, df = df)$value
[1] 0.03297456
2 * (1 - pt(abs(t.test(men, women)$statistic), t.test(men, women)$parameter) )
         t 
0.03297456 

Bayes Factor

Comparaison de deux modèles:

  • \( H_{0}: \mu_{1} = \mu_{2} \rightarrow \delta = 0 \)
  • \( H_{1}: \mu_{1} \neq \mu_{2} \rightarrow \delta \neq 0 \)

\[ \underbrace{\dfrac{p(H_{0}|D)}{p(H_{1}|D)}}_{posterior\ odds} = \underbrace{\dfrac{p(D|H_{0})}{p(D|H_{1})}}_{Bayes\ factor} \times \underbrace{\dfrac{p(H_{0})}{p(H_{1})}}_{prior\ odds} \]

\[ \text{evidence}\ = p(D|H) = \int p(\theta|H) p(D|\theta,H) \text{d}\theta \]

L'évidence en faveur d'un modèle correspond à la probabilité marginale d'un modèle (le dénominateur du théorème de Bayes), c'est à dire à la likelihood moyennée sur le prior… Ce qui fait du Bayes Factor un ratio de likelihood pondéré (par le prior).

Bayes Factor, exemple

On lance une pièce 100 fois et on essaye d'estimer la probabilité \( \theta \) (le biais de la pièce) d'obtenir face. On compare deux modèles qui diffèrent par leur prior sur \( \theta \).

\[ \begin{align} \mathcal{M_{1}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(6, 10) \\ \end{align} \]

\[ \begin{align} \mathcal{M_{2}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(20, 12) \\ \end{align} \]

plot of chunk unnamed-chunk-12

Bayes Factor, exemple


\[ \begin{align} \mathcal{M_{1}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(6, 10) \\ \end{align} \]

\[ \begin{align} \mathcal{M_{2}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(20, 12) \\ \end{align} \]


\[ \text{BF}_{12} = \dfrac{p(D|H_{1})}{p(D|H_{2})} = \dfrac{\int p(\theta|H_{1}) p(D|\theta,H_{1}) \text{d}\theta}{\int p(\theta|H_{2}) p(D|\theta,H_{2}) \text{d}\theta} = \dfrac{\int \mathrm{Binomial}(n, \theta)\mathrm{Beta}(6, 10)\text{d}\theta}{\int \mathrm{Binomial}(n, \theta)\mathrm{Beta}(20, 12) \text{d}\theta} \]

Bayes Factor, exemple

Le Bayes Factor est la nouvelle p-value

Attention à ne pas interpréter le BF comme un posterior odds

Le BF est un facteur qui nous indique de combien nos prior odds doivent changer, au vue des données. Il ne nous dit pas quelle est l'hypothèse la plus probable, sachant les données (sauf si les prior odds sont 1:1).

  • \( H_{0} \): there is no such thing as precognition
  • \( H_{1} \): precognition exists

On fait une expérience et on calcule un \( BF_{10} = 27 \). Quelle est la probabilité a posteriori de H1 ?


\[ \underbrace{\dfrac{p(H_{1}|D)}{p(H_{0}|D)}}_{posterior\ odds} = \underbrace{\dfrac{27}{1}}_{Bayes\ factor} \times \underbrace{\dfrac{1}{1000}}_{prior\ odds} = \dfrac{27}{1000} = 0.027 \]

Bayes Factor - test d'hypothèse nulle

library(BayesFactor)
ttestBF(men, women)
Bayes factor analysis
--------------
[1] Alt., r=0.707 : 1.312345 ±0%

Against denominator:
  Null, mu1-mu2 = 0 
---
Bayes factor type: BFindepSample, JZS

Prior JZS, en référence à Jeffreys (1961), et Zellner & Siow (1980).

JZS prior

data.frame(x = c(-10, 10) ) %>%
    ggplot(aes(x = x) ) +
    stat_function(
        fun = dcauchy,
        args = list(location = 0, scale = sqrt(2) / 2), size = 1.5) +
    theme_bw(base_size = 20)

plot of chunk unnamed-chunk-14

\[ \mathrm{Cauchy} = \dfrac{1}{\pi \gamma \Big[ 1 + \big(\frac{x-x_{0}}{\gamma}\big)^{2}\Big]} \]

Double sens

Dans le cadre bayésien, le terme de prior peut faire référence soit:

  • à la probabilité a priori d'un modèle (par rapport à un autre modèle), c'est à dire \( p(M_{i}) \). Voir ce blogpost.

  • ou aux prédictions du modèle, e.g., \( \beta \sim \mathrm{Normal}(0,10) \). Voir ce blogpost.

plot of chunk unnamed-chunk-15

Comparaison de modèles

Deux problèmes récurrents à éviter en modélisation: l'overffitting et l'underfitting. Comment s'en sortir ?

  • Utiliser des priors pour régulariser, pour contraindre le posterior (i.e., accorder moins de poids à la likelihood)
  • Utiliser des critères d'information (e.g., AIC, WAIC)


\[ R^{2} = \frac{\text{var}(\text{outcome}) - \text{var}(\text{residuals})}{\text{var}(\text{outcome})} = 1 - \frac{\text{var}(\text{residuals})}{\text{var}(\text{outcome})} \]

Overfitting

ppnames <- c("afarensis","africanus","habilis","boisei",
        "rudolfensis","ergaster","sapiens")
brainvolcc <- c(438, 452, 612, 521, 752, 871, 1350)
masskg <- c(37.0, 35.5, 34.5, 41.5, 55.5, 61.0, 53.5)

d <- data.frame(species = ppnames, brain = brainvolcc, mass = masskg)

d %>%
    ggplot(aes(x = mass, y = brain, label = species) ) +
    geom_point() +
    ggrepel::geom_label_repel(hjust = 0, nudge_y = 50, size = 5) +
    theme_bw(base_size = 20) + xlim(30, 70)

plot of chunk unnamed-chunk-16

Overfitting

mod1.1 <- lm(brain ~ mass, data = d)
(var(d$brain) - var(residuals(mod1.1) ) ) / var(d$brain)
[1] 0.490158
mod1.2 <- lm(brain ~ mass + I(mass^2), data = d)
(var(d$brain) - var(residuals(mod1.2) ) ) / var(d$brain)
[1] 0.5359967
mod1.3 <- lm(brain ~ mass + I(mass^2) + I(mass^3), data = d)
(var(d$brain) - var(residuals(mod1.3) ) ) / var(d$brain)
[1] 0.6797736

Overfitting

plot of chunk unnamed-chunk-19

Underfitting

\[ \begin{align} v_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i} &= \alpha \\ \end{align} \]

mod1.7 <- lm(brain ~ 1, data = d)

plot of chunk unnamed-chunk-21

Théorie de l'information

On voudrait mesurer la distance entre notre modèle et le full model (i.e., la “réalité”, la nature)…

De combien notre incertitude est-elle réduite en apprenant un outcome supplémentaire ? Cette réduction est la définition de l'information.

Quelle mesure pour l'incertitude ? S'il existe \( n \) événements possibles, et que chaque événement \( i \) a pour probabilité \( p_{i} \), alors une mesure de l'incertitude est:

\[ H(p) = - \text{E log}(p_{i}) = - \sum_{i=1}^{n}p_{i} \text{log}(p_{i}) \]

En d'autres termes, l'incertitude contenue dans une distribution de probabilités est la log-probabilité moyenne d'un événement.

Incertitude

Exemple de prédiction météorologique. Si on imagine que la probabilité qu'il pleuve ou qu'il fasse beau (à Grenoble) sont \( p_{1}=0.3 \) et \( p_{2}=0.7 \).

Alors, \( H(p) = - (p_{1} \text{log}(p_{1}) + p_{2} \text{log}(p_{2})) \approx 0.61 \).

p <- c(0.3, 0.7)
- sum(p * log(p) )
[1] 0.6108643

Imaginons que nous habitions à Abu Dabi, et que la probabilité qu'il pleuve ou qu'il fasse beau sont \( p_{1}=0.01 \) et \( p_{2}=0.99 \).

p <- c(0.01, 0.99)
- sum(p * log(p) )
[1] 0.05600153

Divergence

On a donc un moyen de quantifier l'incertitude. Comment utiliser cette mesure pour quantifier la distance entre notre modèle et la réalité ?

Divergence: incertitude ajoutée par l'utilisation d'une distribution de probabilités pour décrire… une autre distribution de probabilités (Kullback-Leibler divergence, ou entropie relative).

\[ D_{KL}(p,q) = \sum_{i} p_{i}\big(\text{log}(p_{i}) - \text{log}(q_{i})\big) = \sum_{i} p_{i} \text{log}\bigg(\frac{p_{i}}{q_{i}}\bigg) \]

Divergence

\[ D_{KL}(p,q) = \sum_{i} p_{i}\big(\text{log}(p_{i}) - \text{log}(q_{i})\big) = \sum_{i} p_{i} \text{log}\bigg(\frac{p_{i}}{q_{i}}\bigg) \]

Par exemple, supposons que la véritable distribution de nos événements (soleil vs pluie) soit \( p_{1}=0.3 \) et \( p_{2}=0.7 \). Si nous pensons plutôt que ces événements arrivent avec une probabilité \( q_{1}=0.25 \) et \( q_{2}=0.75 \), quelle quantité d'incertitude avons-nous ajoutée ?

p <- c(0.3, 0.7)
q <- c(0.25, 0.75)

sum(p * log(p / q) )
[1] 0.006401457
sum(q * log(q / p) )
[1] 0.006164264

Entropie croisée et divergence

Entropie croisée: \( H(p,q) = \sum_{i} p_{i} \log (q_{i}) \)

sum(p * (log(q) ) )

La Divergence est définie comme l'entropie additionnelle ajoutée en utilisant \( q \) pour décrire \( p \).

\[ \begin{align} D_{KL}(p,q) &= H(p,q) - H(p) \\ &= - \sum_{i} p_{i} \log(q_{i}) - \big( - \sum_{i} p_{i} \log(p_{i}) \big) \\ &= - \sum_{i} p_{i} \big(\log(q_{i}) - \log(p_{i}) \big) \\ \end{align} \]

- sum (p * (log(q) - log(p) ) )
[1] 0.006401457

Vers la déviance...

OK, mais nous ne connaissons pas la distribution target (la réalité), à quoi cela peut donc nous servir ?

Astuce: si nous comparons deux modèles, \( q \) et \( r \), pour approximer \( p \), nous allons comparer leurs divergences… Et donc \( \text{E} \ \text{log}(p_{i}) \) sera la même quantité pour les deux modèles !

plot of chunk unnamed-chunk-27

Vers la déviance...

On peut donc utiliser \( \text{E} \ \text{log}(q_{i}) \) et \( \text{E} \ \text{log}(r_{i}) \) comme estimateurs de la distance entre chaque modèle et notre distribution cible. On a donc seulement besoin de la log-probabilité moyenne des modèles. Comme on ne connaît pas la distribution cible, cela veut dire qu'on ne peut pas interpréter cette quantité en termes absolus mais seulement en termes relatifs. Ce qui nous intéresse c'est \( \text{E} \ \text{log}(q_{i}) - \text{E} \ \text{log}(r_{i}) \).

plot of chunk unnamed-chunk-28

Déviance

Pour approximer la valeur de \( \text{E} \ \log(q_{i}) \), on peut utiliser la déviance d'un modèle, qui est une mesure du fit relatif du modèle.

\[ D(q) = -2 \sum_{i} \log(q_{i}) \]

où \( i \) indice chaque observation et \( q_{i} \) est la likelihood de chaque observation.

d$mass.s <- scale(d$mass)
mod1.8 <- lm(brain ~ mass.s, data = d)

-2 * logLik(mod1.8) # compute deviance
'log Lik.' 94.92499 (df=3)

Déviance

# extracting model's coefficients

alpha <- coef(mod1.8)[1]
beta <- coef(mod1.8)[2]

# computing the log-likelihood

ll <- sum(dnorm(
    d$brain,
    mean = alpha + beta * d$mass.,
    sd = sd(residuals(mod1.8) ),
    log = TRUE)
    )

# computing the deviance

(-2) * ll
[1] 95.00404

In-sample and out-of-sample

La déviance a le même problème que le \( R^{2} \), lorsqu'elle est calculée sur l'échantillon observé. Dans ce cas, on l'appelle déviance in-sample.

Si on est intéressé par les capacités de prédiction de notre modèle, nous pouvons calculer la déviance du modèle sur de nouvelles données… qu'on appellera dans ce cas déviance out-of-sample. Cela revient à se demander si notre modèle est performant pour prédire de nouvelles données.

Imaginons que nous disposions d'un échantillon de taille \( N \), que nous appellerons échantillon d'apprentissage (training). Nous pouvons calculer la déviance du modèle sur cet échantillon (\( D_{train} \) ou \( D_{in} \)). Si nous acquérons ensuite un nouvel échantillon de taille \( N \) issu du même processus de génération de données (que nous appellerons échantillon de test), nous pouvons calculer une déviance sur ce nouvel échantillon, en utilisant les paramètres estimés avec l'échantillon d'entraînement (que nous appellerons \( D_{test} \) ou \( D_{out} \)).

In sample and out of sample deviance

\[ \begin{align} y_{i} &\sim \mathrm{Normal}(\mu_{i},1) \\ \mu_{i} &= (0.15) x_{1,i} - (0.4) x_{2,i} \\ \end{align} \]

On a réalisé ce processus 10.000 fois pour cinq modèles de régression linéaire de complexité croissante. Les points bleus représentent la déviance calculée sur l'échantillon d'apprentissage et les points noirs la déviance calculée sur l'échantillon de test.

Régularisation

Une autre manière de lutter contre l'overfitting est d'utiliser des priors sceptiques qui vont venir ralentir l'apprentissage réalisé sur les données (i.e., accorder plus de poids au prior).

plot of chunk unnamed-chunk-31

Régularisation

Comment décider de la précision du prior ? Est-ce que le prior est “assez” régularisateur ou pas ?

On peut diviser le jeu de données en deux parties (training et test) afin de choisir le prior qui produit la déviance out-of-sample la plus faible. On appelle cette stratégie la cross-validation.

Critères d'information

On mesure ici la différence entre la déviance in-sample (en bleu) et la déviance out-of-sample (en noir). On remarque que la déviance out-of-sample est presque exactement égale à la déviance in-sample, plus deux fois le nombre de paramètres du modèle…

Akaike information criterion

L'AIC fournit une approximation de la déviance out-of-sample:

\[ \text{AIC} = D_{train} + 2p \]

où \( p \) est le nombre de paramètres libres (i.e., à estimer) dans le modèle. L'AIC donne donc une approximation des capacités de prédiction du modèle.

plot of chunk unnamed-chunk-32

NB: l'AIC fonctionne bien uniquement quand le nombre d'observations \( N \) est largement supérieur au nombre de paramètres \( p \). Dans le cas contraire, on utilise plutôt l'AICc (corrected AIC, voir Burnham & Anderson, 2002; 2004).

Deviance information criterion

Une autre condition de l'AIC est que les priors soient plats (i.e., peu informatifs) ou dépassés par la likelihood (on a beaucoup de données). Le DIC est un indice qui ne requiert pas cette condition, en s'accommodant de priors informatifs.

Le DIC est calculé à partir de la distribution a posteriori de la déviance calculée sur l'échantillon d'apprentissage (i.e., \( D_{train} \)).

\[ \text{DIC} = \bar{D} + (\bar{D} - \widehat{D} ) = \widehat{D} + p_{D} \]

où \( \bar{D} \) est la moyenne de la distribution a posteriori \( D \) calculée pour chaque valeur de paramètre échantillonnée, et \( \widehat{D} \) la déviance calculée à la moyenne de la distribution a posteriori. La différence \( \bar{D} - \widehat{D} = p_{D} \) est analogue au nombre de paramètres utilisé dans le calcul de l'AIC.

La fonction rethinking::DIC permet de le calculer.

Widely applicable information criterion

Une condition d'application de l'AIC et du DIC est que la distribution a posteriori soit une distribution gaussienne multivariée. Le WAIC relache cette condition, en étant souvent plus précis que le DIC.

Un aspect important du WAIC est qu'il est dit pointwise, c'est à dire qu'il considère l'imprécision de prédiction point par point (donnée par donnée), indépendemment pour chaque observation.

On va commencer par calculer la log-pointwise-predictive-density (lppd), définie par:

\[ \text{lppd} = \sum_{i=1}^{N} \log \Pr(y_{i}) \]

où \( \Pr(y_{i}) \) représente la likelihood moyenne de l'observation \( i \) dans le training sample.

En français: la log-densité prédictive point par point est la somme du log de la likelihood moyenne de chaque observation. Il s'agit de l'analogue point par point de la déviance, moyenné sur toute la distribution a posteriori.

Widely applicable information criterion

library(rethinking)
data(cars)

m <- rethinking::map(alist(
        dist ~ dnorm(mu, sigma),
        mu <- a + b * speed,
        a ~ dnorm(0, 100),
        b ~ dnorm(0, 10),
        sigma ~ dunif(0, 30) ),
        data = cars
    )

post <- extract.samples(m, n = 1e3)

head(post)
          a        b    sigma
1 -20.74432 4.068218 16.08980
2 -20.34022 4.146951 12.95898
3 -17.60990 3.990101 17.70716
4 -17.00744 3.963182 17.01450
5 -26.66950 4.280532 16.19768
6 -19.28354 3.969091 15.49696

Widely applicable information criterion

\[ \text{lppd} = \sum_{i=1}^{N} \log \Pr(y_{i}) \]

ll <- sapply(
    1:nrow(post),
    function(s) {

        mu <- post$a[s] + post$b[s] * cars$speed
        dnorm(cars$dist, mu, post$sigma[s], log = TRUE)

        }
    )
lppd <- sapply(1:nrow(cars), function(i) log_sum_exp(ll[i, ]) - log(nrow(post) ) )

sum(lppd)
[1] -206.5098

Widely applicable information criterion

La deuxième partie du calcul du WAIC est le nombre de paramètres effectif, \( p_{WAIC} \). On définit \( V(y_{i}) \) comme la variance de la log-likelihood pour chaque observation \( i \) de l'échantillon d'entraînement.

\[ p_{\text{WAIC}} = \sum_{i=1}^{N} V(y_{i}) \]

pWAIC <- sapply(1:nrow(cars), function(i) var(ll[i, ]) )

Ensuite, le WAIC est défini par:

\[ \text{WAIC} = -2(\text{lppd} - p_{\text{WAIC}}) \]

(WAIC <- (-2) * (sum(lppd) - sum(pWAIC) ) )
[1] 421.0726

Le WAIC est aussi un estimateur de la déviance out-of-sample. La fonction WAIC du package rethinking permet de le calculer pour des modèles construits avec map ou map2stan.

Critères d'information et régularisation

Le DIC et le WAIC peuvent être conceptualisés (au même titre que l'AIC) comme des approximations de la déviance out-of-sample. On remarque que le WAIC produit des approximations plus précises que le DIC, et que l'usage combiné des critères d'information et des priors régularisateurs permet une meilleure approximation de la déviance out-of-sample.

Et après ?

  • Sélection de modèle. On choisit le meilleur modèle un utilisant un des outils présentés et on base nos conclusions sur les paramètres estimés par ce meilleur modèle.

  • Comparaison de modèles. On utilise DIC/WAIC mais aussi les outils de posterior predictive checking discutés précédemment, pour chaque modèle.

  • Moyennage de modèles. On va construire des posterior predictive checks qui exploitent ce qu'on sait des capacités de prédiction de chaque modèle (via le DIC/WAIC).

Comparaison de modèles

Nous allons essayer de prédire les kg par gramme de lait (kcal.per.g) avec les prédicteurs neocortex et le logarithme de mass. Nous allons ensuite fitter 4 modèles qui correspondent aux 4 combinaisons possibles de prédicteurs et les comparer en utilisant le WAIC.

data(milk)

d <- milk[complete.cases(milk), ] # removing NAs
d$neocortex <- d$neocortex.perc / 100 # rescaling explanatory variable

head(d)
              clade            species kcal.per.g perc.fat perc.protein
1     Strepsirrhine     Eulemur fulvus       0.49    16.60        15.42
6  New World Monkey Alouatta seniculus       0.47    21.22        23.58
7  New World Monkey         A palliata       0.56    29.66        23.46
8  New World Monkey       Cebus apella       0.89    53.41        15.80
10 New World Monkey         S sciureus       0.92    50.58        22.33
11 New World Monkey   Cebuella pygmaea       0.80    41.35        20.85
   perc.lactose mass neocortex.perc neocortex
1         67.98 1.95          55.16    0.5516
6         55.20 5.25          64.54    0.6454
7         46.88 5.37          64.54    0.6454
8         30.79 2.51          67.64    0.6764
10        27.09 0.68          68.85    0.6885
11        37.80 0.12          58.85    0.5885

Comparaison de modèles

a.start <- mean(d$kcal.per.g)
sigma.start <- log(sd(d$kcal.per.g) )

mod2.1 <- rethinking::map(alist(
        kcal.per.g ~ dnorm(a, exp(log.sigma) ) ),
        data = d,
        start = list(a = a.start, log.sigma = sigma.start) )

mod2.2 <- rethinking::map(alist(
        kcal.per.g ~ dnorm(mu, exp(log.sigma) ),
        mu <- a + bn * neocortex),
        data = d,
        start = list(a = a.start, bn = 0, log.sigma = sigma.start) )

mod2.3 <- rethinking::map(alist(
        kcal.per.g ~ dnorm(mu, exp(log.sigma) ),
        mu <- a + bm * log(mass) ),
        data = d,
        start = list(a = a.start, bm = 0, log.sigma = sigma.start) )

mod2.4 <- rethinking::map(alist(
        kcal.per.g ~ dnorm(mu, exp(log.sigma) ),
        mu <- a + bn * neocortex + bm * log(mass) ),
        data = d,
        start = list(a = a.start, bn = 0, bm = 0, log.sigma = sigma.start) )

Comparaison de modèles

(milk.models <- compare(mod2.1, mod2.2, mod2.3, mod2.4) )
        WAIC pWAIC dWAIC weight   SE  dSE
mod2.4 -14.7   5.0   0.0   0.92 7.91   NA
mod2.1  -8.4   1.8   6.3   0.04 4.52 7.41
mod2.3  -7.8   3.0   6.8   0.03 5.67 5.35
mod2.2  -6.2   2.9   8.4   0.01 4.28 7.79
plot(milk.models, SE = TRUE, dSE = TRUE)

plot of chunk unnamed-chunk-40

Akaike's weights

Le poids d'un modèle est une estimation de la probabilité que ce modèle fera les meilleures prédictions possibles sur un nouveau jeu de données, conditionnellement au set de modèles considéré.

\[ w_{i} = \frac{\exp (-\frac{1}{2}\text{dWAIC}_{i})}{\sum_{j=1}^{m} \exp (-\frac{1}{2}\text{dWAIC}_{j})} \]

Cette formule “remet” simplement le WAIC à l'échelle des probabilités. Le modèle mod2.1 a un poids de \( 0.92 \), qui le place en tête du jeu de modèles. N'oublions cependant pas que nous disposons seulement de 12 observations…

Comparaison des paramètres estimés par nos modèles

coeftab(mod2.1, mod2.2, mod2.3, mod2.4)
          mod2.1  mod2.2  mod2.3  mod2.4 
a            0.66    0.35    0.71   -1.09
log.sigma   -1.79   -1.80   -1.85   -2.16
bn             NA    0.45      NA    2.79
bm             NA      NA   -0.03   -0.10
nobs           17      17      17      17
plot(coeftab(mod2.1, mod2.2, mod2.3, mod2.4) )

plot of chunk unnamed-chunk-41

Moyennage de modèles

Pourquoi ne conserver uniquement le premier modèle et oublier les autres ? Une autre stratégie consisterait à pondérer les prédictions des modèles par leurs poids respectifs. C'est ce qu'on appelle le moyennage de modèles (model averaging).

  • Calculer le WAIC de chaque modèle
  • Calculer le poids de chaque modèle
  • Simuler des données à partir de chaque modèle
  • Combiner ces valeurs simulées dans un ensemble de prédictions, en utilisant les poids des modèles comme proportions
nc.seq <- seq(from = 0.5, to = 0.8, length.out = 30)

d.predict <- list(
    kcal.per.g = rep(0, 30), # empty outcome
    neocortex = nc.seq,      # sequence of neocortex
    mass = rep(4.5, 30)      # average mass
    )

milk.ensemble <- ensemble(mod2.1, mod2.2, mod2.3, mod2.4, data = d.predict, refresh = 0)
mu <- apply(milk.ensemble$link, 2, mean)
mu.HPDI <- t(apply(milk.ensemble$link, 2, HPDI) )

Moyennage de modèles

Voici les prédictions de tous les modèles considérés, pondérés par leur poids respectif. Comme le modèle mod2.4 concentrait quasiment tout le poids, il fait sens que cette prédiction moyennée soit similaire aux prédictions du modèle mod2.4.

plot of chunk unnamed-chunk-43

Travaux pratiques

data(Howell1)

d <-
    Howell1 %>%
    mutate(age = scale(age) )

set.seed(1000)
i <- sample(1:nrow(d), size = nrow(d) / 2)

d1 <- d[i, ] # training sample
d2 <- d[-i, ] # test sample

Nous avons maintenant deux dataframes, de 272 lignes chacune. On va utiliser d1 pour fitter nos modèles et d2 pour les évaluer.

Travaux pratiques

Soit \( h_{i} \) les valeurs de taille et \( x_{i} \) les valeurs centrées d'âge, sur la ligne \( i \). Construisez les modèles suivants avec d1, en utilisant rethinking::map et des priors faiblement régularisateurs.

\[ \begin{align} \mathcal{M_{1}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} \\ \mathcal{M_{2}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} \\ \mathcal{M_{3}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} + \beta_{3} x_{i}^{3} \\ \mathcal{M_{4}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} + \beta_{3} x_{i}^{3} + \beta_{4} x_{i}^{4} \\ \mathcal{M_{5}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} + \beta_{3} x_{i}^{3} + \beta_{4} x_{i}^{4} + \beta_{5} x_{i}^{5} \\ \mathcal{M_{6}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} + \beta_{3} x_{i}^{3} + \beta_{4} x_{i}^{4} + \beta_{5} x_{i}^{5} + \beta_{6} x_{i}^{6} \\ \end{align} \]

Travaux pratiques

  1. Comparer ces modèles en utilisant le WAIC. Comparer les rangs des modèles et leurs poids.
  2. Pour chaque modèle, produire un plot de la moyenne estimée et l'intervalle de confiance à 97% de la moyenne, surimposée aux données brutes. Comment ces prédictions diffèrent-elles selon les modèles ?
  3. Faire un plot des prédictions moyennées sur tous les modèles (sur les trois meilleurs). En quoi ces prédictions diffèrent-elles des prédictions du modèle avec le plus petit WAIC ?
  4. Calculer la déviance out-of-sample pour chaque modèle. Vous pouvez calculer la log-likelihood des tailles en utilisant sum(dnorm(d2$height, mu, sigma, log = TRUE))mu est un vecteur de moyennes prédites, et sigma est l'écart type produit par map.
  5. Comparer les déviances obtenues à la question précédente aux valeurs de WAIC. Basé sur les déviances obtenues, quel modèle fait les meilleures prédictions ? Est-ce que le WAIC est un bon estimateur la déviance out-of-sample ?

Solution - question 1

alpha.start <- mean(d1$height)
sigma.start <- sd(d1$height)

mod3.1 <- rethinking::map(alist(
        height ~ dnorm(mean = mu, sd = sigma),
        mu <- alpha + beta.1*age,
        c(alpha, beta.1) ~ dnorm(mean = 0, sd = 100),
        sigma ~ dunif(min = 0, max = 50) ),
        data = d1,
        start = list(alpha = alpha.start) )

mod3.2 <- rethinking::map(alist(
        height ~ dnorm(mean = mu, sd = sigma),
        mu <- alpha + beta.1*age + beta.2*age^2,
        c(alpha, beta.1, beta.2) ~ dnorm(mean = 0, sd = 100),
        sigma ~ dunif(min = 0, max = 50) ),
        data = d1,
        start = list(alpha = alpha.start) )

Solution - question 1

mod3.3 <- rethinking::map(alist(
        height ~ dnorm(mean = mu, sd = sigma),
        mu <- alpha + beta.1*age + beta.2*age^2 + beta.3*age^3,
        c(alpha, beta.1, beta.2, beta.3) ~ dnorm(mean = 0, sd = 100),
        sigma ~ dunif(min = 0, max = 50) ),
        data = d1,
        start = list(alpha = alpha.start) )

mod3.4 <- rethinking::map(alist(
        height ~ dnorm(mean = mu, sd = sigma),
        mu <- alpha + beta.1*age + beta.2*age^2 + beta.3*age^3 + beta.4*age^4,
        c(alpha, beta.1, beta.2, beta.3, beta.4) ~ dnorm(mean = 0, sd = 100),
        sigma ~ dunif(min = 0, max = 50) ),
        data = d1,
        start = list(alpha = alpha.start) )

Solution - question 1

mod3.5 <- rethinking::map(alist(
        height ~ dnorm(mean = mu, sd = sigma),
        mu <- alpha + beta.1*age + beta.2*age^2 + beta.3*age^3 + 
                beta.4*age^4 + beta.5*age^5,
        c(alpha, beta.1, beta.2, beta.3, beta.4, beta.5) ~ dnorm(mean = 0, sd = 100),
        sigma ~ dunif(min = 0, max = 50) ),
        data = d1,
        start = list(alpha = alpha.start) )

mod3.6 <- rethinking::map(alist(
        height ~ dnorm(mean = mu, sd = sigma),
        mu <- alpha + beta.1*age + beta.2*age^2 + beta.3*age^3 +
                beta.4*age^4 +
                beta.5*age^5 + beta.6*age^6,
        c(alpha, beta.1, beta.2, beta.3, beta.4, beta.5, beta.6) ~dnorm(mean = 0, sd = 100),
        sigma ~ dunif(min = 0, max = 50) ),
        data = d1,
        start = list(alpha = alpha.start) )

Solution - question 1

compare(mod3.1, mod3.2, mod3.3, mod3.4, mod3.5, mod3.6)
         WAIC pWAIC dWAIC weight    SE   dSE
mod3.4 1926.0   5.6   0.0   0.57 25.44    NA
mod3.5 1927.5   6.4   1.4   0.28 25.32  1.15
mod3.6 1928.6   7.8   2.6   0.16 25.04  2.89
mod3.3 1952.3   5.4  26.2   0.00 24.19 10.85
mod3.2 2150.0   5.3 224.0   0.00 22.61 26.72
mod3.1 2395.4   3.4 469.4   0.00 22.94 31.09

Solution - question 2

n.trials <- 1e4
age.seq <- seq(from = -2, to = 3.5, length.out = 58)
prediction.data <- data.frame(age = age.seq)

computeMu <- function(model, data, n.trials) {

    mu <- link(fit = model, data = data, n = n.trials)
    return(mu)

}

computeMuMean <- function(mu) {

    mu.mean <- apply(X = mu, MARGIN = 2, FUN = mean)
    return(mu.mean)

}

computeMuHPDI <- function(mu) {

    mu.HPDI <- apply(X = mu, MARGIN = 2, FUN = HPDI, prob = 0.97)
    return(mu.HPDI)

}

Solution - question 2

simulateHeights <- function(model, prediction.data) {

    simulated.heights <- sim(fit = model, data = prediction.data)
    return(simulated.heights)

}

plotResults <- function(model, prediction.data, original.data, n.trials) {

    mu <- computeMu(model, prediction.data, n.trials)
    mu.mean <- computeMuMean(mu)
    mu.HPDI <- computeMuHPDI(mu)
    simulated.heights <- simulateHeights(model = model, prediction.data = prediction.data)
    simulated.heights.HPDI <- apply(X = simulated.heights, MARGIN = 2, FUN = HPDI)
    plot(height ~ age, data = original.data, col = "steelblue", pch = 16)
    lines(x = prediction.data$age, y = mu.mean, lty = 2)
    lines(x = prediction.data$age, y = mu.HPDI[1, ], lty = 2)
    lines(x = prediction.data$age, y = mu.HPDI[2, ], lty = 2)
    shade(object = simulated.heights.HPDI, lim = prediction.data$age)

}

Solution - question 2

plotResults(
    model = mod3.1, prediction.data = prediction.data,
    original.data = d1, n.trials = n.trials
    )

plotResults(
    model = mod3.2, prediction.data = prediction.data,
    original.data = d1, n.trials = n.trials
    )

plotResults(
    model = mod3.3, prediction.data = prediction.data,
    original.data = d1, n.trials = n.trials
    )

plotResults(
    model = mod3.4, prediction.data = prediction.data,
    original.data = d1, n.trials = n.trials
    )

plotResults(
    model = mod3.5, prediction.data = prediction.data,
    original.data = d1, n.trials = n.trials
    )

plotResults(
    model = mod3.6, prediction.data = prediction.data,
    original.data = d1, n.trials = n.trials
    )

Solution - question 3

h.ensemble <- ensemble(mod3.4, mod3.5, mod3.6, data = list(age = age.seq) )
mu.mean <- apply(h.ensemble$link, 2, mean)
mu.ci <- t(apply(h.ensemble$link, 2, HPDI) )
height.ci <- t(apply(h.ensemble$sim, 2, HPDI) )

ggplot(data = d1, aes(x = as.numeric(age), y = height) ) +
    geom_point(size = 2) +
    geom_line(data = data.frame(age = age.seq, height = mu.mean) ) +
    geom_ribbon(
        data = data.frame(mu.ci), inherit.aes = FALSE,
        aes(x = age.seq, ymin = mu.ci[, 1], ymax = mu.ci[, 2]), alpha = 0.2) +
    geom_ribbon(
        data = data.frame(height.ci), inherit.aes = FALSE,
        aes(x = age.seq, ymin = height.ci[, 1], ymax = height.ci[, 2]), alpha = 0.1) +
    theme_bw(base_size = 20) + xlab("age") + xlim(-2, 3) + ylim(60, 190)

plot of chunk unnamed-chunk-52

Solution - question 4

# model 1
coefs <- coef(mod3.1)
mu <- coefs["alpha"] + coefs["beta.1"] * d2$age
log.likelihood <- sum(dnorm(x = d2$height, mean = mu, sd = coefs["sigma"], log = TRUE) )
dev.mod3.1 <- -2 * log.likelihood

# model 2
coefs <- coef(mod3.2)
mu <- coefs["alpha"] + coefs["beta.1"] * d2$age + coefs["beta.2"] * (d2$age)^2
log.likelihood <- sum(dnorm(x = d2$height, mean = mu, sd = coefs["sigma"], log = TRUE) )
dev.mod3.2 <- -2 * log.likelihood

# model 3
coefs <- coef(mod3.3)
mu <- coefs["alpha"] + coefs["beta.1"] * d2$age + coefs["beta.2"] * (d2$age)^2 + coefs["beta.3"] * (d2$age)^3
log.likelihood <- sum(dnorm(x = d2$height, mean = mu, sd = coefs["sigma"], log = TRUE) )
dev.mod3.3 <- -2 * log.likelihood

Solution - question 4

# model 4
coefs <- coef(mod3.4)
mu <- coefs["alpha"] + coefs["beta.1"]*d2$age + coefs["beta.2"]*(d2$age)^2 + coefs["beta.3"]*(d2$age)^3 + coefs["beta.4"]*(d2$age)^4
log.likelihood <- sum(dnorm(x = d2$height, mean = mu, sd = coefs["sigma"], log = TRUE) )
dev.mod3.4 <- -2 * log.likelihood

# model 5
coefs <- coef(mod3.5)
mu <- coefs["alpha"] + coefs["beta.1"]*d2$age + coefs["beta.2"]*(d2$age)^2 + coefs["beta.3"]*(d2$age)^3 + coefs["beta.4"]*(d2$age)^4 + coefs["beta.5"]*(d2$age)^5
log.likelihood <- sum(dnorm(x = d2$height, mean = mu, sd = coefs["sigma"], log = TRUE) )
dev.mod3.5 <- -2 * log.likelihood

# model 6
coefs <- coef(mod3.6)
mu <- coefs["alpha"] + coefs["beta.1"]*d2$age + coefs["beta.2"]*(d2$age)^2 + coefs["beta.3"]*(d2$age)^3 + coefs["beta.4"]*(d2$age)^4 + coefs["beta.5"]*(d2$age)^5 + coefs["beta.6"]*(d2$age)^6
log.likelihood <- sum(dnorm(x = d2$height, mean = mu, sd = coefs["sigma"], log = TRUE) )
dev.mod3.6 <- -2 * log.likelihood

Solution - question 5

deviances <- c(dev.mod3.1,dev.mod3.2,dev.mod3.3,dev.mod3.4,dev.mod3.5,dev.mod3.6)
comparison <- compare(mod3.1,mod3.2,mod3.3,mod3.4,mod3.5,mod3.6)
comparison <- as.data.frame(comparison@output)
comparison <- comparison[order(rownames(comparison) ), ]
waics <- comparison$WAIC

data.frame(deviance = deviances, waic = waics) %>%
    gather(type, value) %>%
    mutate(x = rep(1:6, 2) ) %>%
    ggplot(aes(x = x, y = value, colour = type) ) +
    scale_colour_grey() +
    geom_point(size = 2) +
    scale_x_continuous(breaks = 1:6) +
    theme_bw(base_size = 20) + xlab("model") + ylab("Déviance/WAIC")

plot of chunk unnamed-chunk-55